Anomalous Heat Conduction in a Di-atomic One-Dimensional Ideal Gas 
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We provide firm convincing evidence that the energy trans- 
port in a one- dimensional gas of elastically colliding free par- 
ticles of unequal masses is anomalous, i.e. the Fourier Law 
does not hold. Our conclusions are based on the analysis of 
the dependence of the heat current on the number of particles, 
of the internal temperature profile and on the Green-Kubo 
formalism. 

PACS number: 05.45.-a 



Understanding the dynamical origin for the validity 
of the Fourier law of heat conduction in deterministic 
one-dimensional particle chains is one of the oldest and 
most frustrating problems in non-equilibrium statistical 
physics Due to some very basic unresolved issues 

the problem has been a source of many recent publica- 
tions §-§. 

In the absence of analytical results, these papers are 
mainly oriented towards a numerical analysis of the prob- 
lem. However, due to the delicate nature of the questions 
under discussion, numerical results sometimes lead to dif- 
ferent conclusions. This is the case, for example, of the 
Id hard-point particles with alternating masses for which 
opposite conclusions have been reached |||7|,|| . This dis- 
agreement is not extremely surprising since this system 
lies in the foggy region which separates clear, regular 
integrable systems, from the totally chaotic, determin- 
istic, motion. Indeed this system has zero Lyapounov 
exponent and therefore it lacks of the exponential local 
instability which characterizes chaotic systems. On the 
other hand it has been shown |l^] that such systems can 
exhibit Gaussian diffusive behaviour and, more recently 
JnJ, an example has been shown of a system with zero 
Lyapounov exponent which however obeys the Fourier 
law. From the point of view of a general theoretical un- 
derstanding, the fact that the alternating mass problem 
lies in this critical region renders particularly important 
to establish, beyond any reasonable doubt, its conducting 
properties. This is what we set up to do in the present 
paper. 

We consider a one-dimensional gas of interacting par- 
ticles with coordinates q n and momenta p n for which the 
hamiltonian can be written in the form 
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The energy current from site n to site n + 1 is defined 
as j n = {h n+ i,h n } and satisfies the continuity equation 
{d/dt)h n = {H, h n } = j n — j n -i- In particular we focus 
our attention on the ideal gas model of elastically collid- 
ing particles, V(q > 0) = 0, V(q < 0) — oo, with alter- 
nating masses, m<x n -\ = m\ = y/r, rri2n = tcii = 
where the ratio r = TO1/7712 serves as a model parameter. 
We have mainly considered the value r = (v5 — 1 )/2; 
however all the reported numerical results have been 
checked also for several other values of r (0 < r < 1) 
where we found no qualitative distinction. 

We place our system of N particles between two 
stochastic maxwellian heat reservoirs at temperatures Tl 
and Tr (sec || for a description of the reservoir model). 
We chose the temperatures of the reservoirs, Tl = I 
and Tr, = 2, and measure the long-time averaged heat 
current (J) — Umt^. 0O (l/t) f dt'J(t') versus the system 

size N, where J = (1/N) ^2n=i Here we want to 
note that strict equivalence between our definition of the 
heat current (which simply accounts for the energy trans- 
ferred during collisions), and the 'free particle' current 
j n = m n Vn/2 used by some authors, e.g S (which does 
not 'feel' the collisions), is nontrivial in general ||. 

In order to ensure that our results are not affected by 
finite size effects we have put particular care in using 
an efficient numerical scheme which allows to reach high 
N values. Our algorithm, developed for the first time 
in Ref. ||, searches in a partially ordered tree (heap) of 
pre-computed candidates pairs for the next collision and, 
due to this, its requires only log 2 N computer operations 
per collision. As a consequence, we were able to simu- 
late very long chains and we have obtained reliably con- 
verged results for lattices with sizes TV up to 30000. Con- 
vergence has been controlled by checking the constancy 
of the finite-time-averaged heat current (1/t) J Q * dt'J(t r ), 
and to this end simulations for the largest system sizes 
had to be carried on up to 10 12 pair collisions. It is also 
clear that convergence problems suggest to keep far away 
form the r values too close to one or to zero. The analysis 
made in indicated that the range 0.15 < r < 0.6 was 
the most effective in attenuating solitary pulses and the 
value r = 0.2 was chosen. In the present paper we take 
the somehow 'standard' choice r = (v5— l)/2- 

Validity of Fourier law implies the scaling J oc VT cx 
A^ 1 . Our numerical results shown in fig. [l] clearly 
demonstrate instead a different power-law behaviour 
namely J oc A^~ a with a sw 0.745 ± 0.005 over a very 
large range in N. We have also found that the scaling 
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exponent a does not change appreciably with the mass 
ratio r. For example for ten times smaller value of r 
the asymptotic scaling only sets in later (i.e. for larger 
values of N, see fig. Q). The possibility of a slow con- 
vergence to the asymptotic value might be at the origin 
of the slightly different numerical values for a found in 
previous numerical experiments (a w 0.65 in Hatano [^J, 
a pa 0.83 in Dhar 0). Since the model under considera- 
tion is energy scaling we do not expect any dependence 
of the exponent a on the reservoirs temperatures. 



o.i 



A 
V 



0.01 



1=0.618 t- 

1=0.0618 — x- 
slope=-0.745 



100 



1000 



10000 



FIG. 1. Time- averaged energy current of a system of 
N particles between heat baths at temperatures TL = 1 
and Tr = 2 vs. the size N, at two different mass ratios 
r = mi/m2 =■ (v5— 1)/2 and r — (v5— 1)/20. The suggested 
scaling (J) oc N~ a with a = 0.745 is shown for comparison. 

The above results therefore solve the existing contro- 
versy and clearly show that the alternating mass, Id 
hard point gas does not obey Fourier heat law. We turn 
now to the analysis of other quantities which, besides 
providing additional confirmation of the above conclu- 
sions, illuminate interesting aspects of the heat conduc- 
tion problem. A quantity of main interest is the inter- 
nal local temperature profile T n = (p^/(2m„)) in the 
non-equilibrium steady state for the system placed in 
between the heat reservoirs. First we notice that the 
temperature profile in discrete index variable n is differ- 
ent than the temperature profile in position variable q n 
[Q] since the inverse density dq/dn = (q n +i — qn) is non- 
uniform in non-equilibrium, in fact it is simply propor- 
tional to the temperature due to constancy of pressure 
0. Now, in case of a Fourier law, the thermal conduc- 
tivity k scales with temperature like k cx \[T . Therefore, 
from \/T(dT/dn)dn/dq — const one obtains the temper- 



ature profile T* 1 n = (Tl 
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sive numerical simulations showed that the temperature 
profile in our model converges, for sufficiently large N, to 
a well-defined scaling function T^ cal = f(n/N) which is 
slightly, but significantly different, from the kinetic tem- 
perature profile T*™. This is another evidence for the 
anomalous heat transport and for the non- validity of the 
Fourier law in our system. It should be remarked that 



the convergence to the temperature profile predicted by 
kinetic theory observed in j7|, which has indeed been con- 
sidered as surprising by the author himself, actually does 
not take place. 

A standard theoretical analysis of transport laws is 
based on Kubo formulae (lj,[ll|. However, applicabil- 
ity of Kubo formula in momentum conserving cases, i.e. 
for translationally invariant systems like model (Q), is 
not very clear. This is particularly critical in view of a 
recent claim Q that Kubo formula diverges for a mo- 
mentum conserving lattice with non-vanishing pressure. 
For this latter type of models we have an additional dif- 
ficulty in applying the Kubo formalism, namely, as we 
show below, the result depends not only on the temper- 
ature gradient, but also on other thermodynamic prop- 
erties of the initial non-equilibrium state - i.e. the iso- 
baric case (constant pressure profile) or the isodense case 
(constant density profile). There is no a-priori argument 
which favours either of these two options: the choice de- 
pends on the specific physical situation of interest. For 
instance the steady-state heat current simulation consid- 
ered above (figs. 1,2) is clearly described by the isobaric 
state. Since we want to consider both situations we need 
to revise the derivation of Kubo formula by following the 
time evolution of a general non-equilibrium initial state 
in an isolated system with periodic boundary conditions 
<7jv = qo + N,pn = po- To this end, we prepare the initial 
state in a local-equilibrium state described by the inverse 
temperature profile /3„ and by an additional thermody- 
namic potential 7„ 
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FIG. 2. Temperature profile for T L = 1, Tr = 2, and dif- 
ferent sizes N from 800 to 12800 (dotted-dashed-solid curves), 
as compared to the Fourier law prediction (chain curve). 

This (small) additional term is necessary in order to equi- 
librate the pressure in the isobaric case. Notice that j n 
is undetermined up to an arbitrary additive constant due 
to a gauge invariance of the second term in (^|). In or- 
der to determine the gradient of the 7-potential which 
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is necessary to keep the physical pressure 4> constant (n- 
independent) , we compute the generalized pressure <j>i 

d 

Ml = -^ln^(a)| Q=0 , (3) 
Z[(a) — j' e "£ n C» y ('»+ 1 "«»+" i '™)^»('»+ 1 -«» + " ,i '"))(Ig 

By a simple trick, a shift of one variable qi — > qi + a in 
the integral Zi(a), we find Zi(a) = Zi-\[a) and therefore 

Mi = $_i<^_i = const. (4) 

Writing the physical pressure (force) as </> = — {V'(q n -\-i — 
<?«))ncq = 4>n + In/ Pn, multiplying by f3 n , and taking the 
first difference V/ n := f n — / n _i we obtain the required 
'gradient' 

V 7 „ = 0V/3„. (5) 

In the following we consider two specific cases: (i) The 
initial isodense state with (q n +i — q n )ncq — const. This 
is obtained by putting 7„ = 0. (ii) The initial isobaric 
state with uniform pressure profile. This is obtained by 
specifying the 7-potential according to eq. (||). We note 
again that the isobaric state (ii) is the one which is rel- 
evant for the steady non-equilibrium state of a system 
in contact with heat reservoirs. Carefully repeating the 
first few steps in the derivation of the Green-Kubo for- 
mula (following Ref. Jl3|]) we arrive at the very general 
linear response formula 

(A(t) - A(t )) acq = f dt'(A(t') ^(V(3 n j n + V7„W„))eq 

where v n — q n are the particles' velocities. In the 
last step we have assumed that we are close to equi- 
librium (V/3„ and Vj n small) so that the RHS can be 
evaluated in the corresponding equilibrium state () oq . 
Let us now consider the periodic temperature profile 
f3 n = [3 + e sin(2irkn/N), and compute the total heat 
that has been transported in time t from warm to cold 
regions of the lattice, namely Q(t) — f*dt'Jk{t) where 

Jfc := iV -1 ^2^=0 cos (2tt kn/N)j n . Inserting Q for A and 
using (||) (case (i) is obtained by formally setting <fi = 0) 
we obtain 

(Q(t))„cq = ^ / dt'(t - t')(J k (t){J k + (f>V k )) cq . (6) 

^ v Jo 

where V k := N^ 1 J2n=o cos(2Trkn/N)v n . We see that 
the growth of the transported heat (Q(t)) ncq is given by 
the generalized correlation function C k (t) — Cjjit) + 
4>Cjv{t), where Cjj{t) = (Jfc(t)J fe ) eq and C JV (t) = 
(Jk(t)Vk)eq- In the isodense case (i) expression (g) re- 
duces to the usual current autocorrelation function only. 

In the case of Fourier law, we expect initial linear 
growth (Q(t)) noq ~ cut, while for the ballistic heat trans- 
port we expect quadratic growth (Q(t)} ncq oc t 2 (this has 



been confirmed numerically for the integrable gas of equal 
masses r = I). However, in a generic system with mo- 
mentum conservation and non- vanishing pressure <p 7^ 0, 
like our diatomic hard-point gas, we find qualitatively dif- 
ferent behavior in cases (i) and (ii). For example, due to 
the result J(|, Cjj(t) has a plateau oc 4> 2 and the trans- 
port is ballistic in the isodense case, while in the iso- 
baric case the second term, Cjy(i), compensates for the 
plateau and yields a much slower increase of the trans- 
ported heat. In this latter case, independent numerical 
computations of (Q{t)) ncq and of Ck(t) for N up to 32768 
shown in fig. ||give (Q(t)} ncq oc t v with v w 1.255, which 
is still clearly super-diffusive, and directly validate the 
formula (||). 
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FIG. 3. Transported heat Q(t) in an isolated system of 
size iV = 2048 obtained by starting from a non-equilibrium 
isobaric initial state (circles) with e = 0.2. For comparison we 
show the corresponding equilibrium averaged Kubo-like ex- 
pressions (§) for N = 2048 (dashed) and for N = 32768 (solid 
curve, multiplied by 16 to account for scaling (Q(t)) oc 1/N). 
The dashed-dotted line has a slope 1.255 and gives the best fit 
in the range 20 < t < 2000. The corresponding data for the 
isodense case are shown in the inset with the ballistic slope 2. 

In fig. ^ we show the generalized correlation function 
Ck{t) for k = and k = 1 separately. Note that the 
results for k = 1 exhibits some oscillations for longer 
times due to finite size effects, namely due to periodicity 
of the lattice, while k = gives the spatially homoge- 
neous correlation function which has the same long-time 
behavior with weaker finite size effects [however, the case 
k = is not strictly related to Kubo formula @j. We see 
that in the time range where Co(t) and C\(t) match, the 
asymptotic decay is compatible with Ck(t) oc with 
the exponent \i = 2 — v = 0.745 consistent with eq. (|^), 
and satisfying /1 = a. 

This results can be interpreted in the following way. 
In the isodense initial state the initial temperature gra- 
dients drive the heat ballistically in terms of sound waves 
Q. On the other hand, in the isobaric initial state, we 
have density gradients which drive the heat in the op- 
posite direction and almost exactly compensate for the 
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effect of temperature gradients so that the net effect is a 
sub-ballistic, but still super-diffusive, energy transport. 
In order to illustrate the mechanism of ballistic energy 
transport we show in fig. ^| the spatio-temporal current- 
current correlation function Cjj(n,t) — (joj n {t)) C q which 
exhibits clear ballistic tongues along the lines n — ±c s i 
where c s = 1.78. 
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FIG. 5. The spatio-temporal current-current correlation 
function Cjj(n,t) at temperature 1//3 = 1 on a lattice of size 
N — 1024 is shown with twenty different shades of greyness 
spaced equidistantly from 10 -4 to 4.0 in logarithmic scale. 
The zig-zag solid line indicates the peak ballistic sound-wave 
contribution moving with a uniform sound velocity c s — 1.78. 



FIG. 4. The generalized time correlation function (see 
text) computed with canonical average at j3 = 1 for two sys- 
tem sizes N = 1024 and N = 8192 and for the zeroth k = 
(dashed/solid curve) and the first k = 1 (thin curve) spatial 
Fourier mode. Note the t 
range 20 < t < 200 (for N 
we see finite size effects (e.g. we have periodic oscillations for 
k — 1 due to transversals of sound waves) . 



745 decay (dashed line) in the 
= 8192) whereas for longer times 



In this paper we have discussed the thermal conduct- 
ing properties of a one dimensional hard point gas with 
alternating masses. This problem has a long history and 
recently several numerical results have appeared which 
lead to contrasting conclusions on whether Fourier law is 
obeyed or not. In the latter case, different values have 
been found for the rate of divergence of coefficient of ther- 
mal conductivity as a function of the particles number. 
Indeed, the slow convergence properties which sometimes 
characterize this problem suggest particular care in the 
interpretation of numerical findings. 

Here we have presented accurate numerical analysis, 
made possible by a powerful integration scheme, which 
allows us to establish definite convincing evidence that 
the system does not obey the Fourier law. Moreover, by 
considering a typical mass ratio r = (v5 — l)/2, we have 
found that the asymptotic scalings: (i) steady-state heat 
current between heat baths (J) oc N~ a , (ii) heat trans- 
ported within a non-equillibrium isobaric initial state 
in isolated system (Q(t)) oc t 2 ~ a , and (iii) generalized 
current- velocity correlation in equilibrium state Ck{t) oc 
t~ a , are described by just one exponent a = 0.745. 
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